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ABSTRACT 

We demonstrate a new Bayesian technique to invert color-magnitude dia- 
grams of main sequence and white dwarf stars to reveal the underlying cluster 
properties of age, distance, metallicity, and line-of-sight absorption, as well as 
individual stellar masses. The advantages our technique has over traditional 
analyses of color-magnitude diagrams are objectivity, precision, and explicit de- 
pendence on prior knowledge of cluster parameters. Within the confines of a 
given set of often-used models of stellar evolution, the initial-final mass relation, 
and white dwarf cooling, and assuming photometric errors that one could reason- 
ably achieve with the Hubble Space Telescope, our technique yields exceptional 
precision for even modest numbers of cluster stars. For clusters with 50 to 400 
members and one to a few dozen white dwarfs, we find typical internal errors of 
(T([Fe/H]) < 0.03 dex, a(m-Mv) < 0.02 mag, and a(Av) < 0.01 mag. We derive 
cluster white dwarf ages with internal errors of typically only 10% for clusters 
with only three white dwarfs and almost always < 5% with ten white dwarfs. 
These exceptional precisions will allow us to test white dwarf cooling models and 
standard stellar evolution models through observations of white dwarfs in open 
and globular clusters. 

Subject headings: open clusters and associations: general — stars: evolution — 
white dwarfs 
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1. Introduction 

White dwarf cooling theory currently provides the most reliable age for the Galactic disk 
(Winget et al. 1987; Oswalt et al. 1996; Leggett, Ruiz, & Bergeron 1998; Knox, Hawkins, 
& Hambly 1999), whereas main sequence stellar evolution provides the most reliable age 
for the Galactic halo (e.g., Salaris & Weiss 2002; Krauss & Ghaboyer 2003). In order to 
understand the detailed formation sequence of the Galactic components, as well as the local 
satellite galaxies, these two time scales need to be placed on the same absolute age system. 
The only current empirical approach available to inter-calibrate these two age systems is 
to derive white dwarf (WD) cooling ages and main sequence turn-off (MSTO) ages for a 
number of Galactic star clusters over a wide range of ages and metallicities. Much of the 
WD age dating work has been necessarily limited to nearby open clusters (Claver 1995; von 
Hippel, Gilmore, & Jones 1995; Richer et al. 1998; von Hippel & Gilmore 2000; Claver et al. 
2001; von Hippel 2005, hereafter paper 1) that are young or of intermediate age, since old 
WDs are faint. Hansen et al. (2002) extended WD age studies to one globular cluster (NGG 
6121 = M4). They derived a precise WD age, but with large systematic uncertainties due 
to as-yet uncalibrated physical effects in the coolest WDs (Fontaine, Brassard, & Bergeron 
2001). 

Even though the Hubble Space Telescope may be nearing the end of its hfetime, it has 
made collecting these deep observations of WDs in open and globular clusters possible. At 
least two more open clusters (NGG 2360, NGG 2660) and one more globular cluster (NGG 
6397) have been observed with HST to sufficient depth, and those results will be forthcoming. 
The large number of 8-lOm telescopes now available make it possible to observe a few more 
open clusters to sufficient depth for the WD technique, and the next decade should see 
20-30m telescopes, which will make these studies substantially easier. 

While the instrumentation has been improving and there has been steady work on 
improving WD cooling and traditional main sequence stellar evolutionary models, there have 
not been sufficient advances in the statistical machinery available to compare star cluster 
observations with those models, particularly for WDs. In this paper we present the first phase 
of our effort to develop this statistical machinery. Specifically, we present a new Bayesian 
technique that has the ability to objectively incorporate all our prior knowledge, including 
stellar evolution, star cluster properties, and data quality estimates, while comparing data 
for each cluster to any available theoretical model. We chose to employ a Bayesian approach 
precisely because so much is known about stellar evolution and star clusters, and because 
this approach allows us to test how cluster properties depend on the input models or model 
ingredients. 

The power of the Bayesian approach is impressive, and we show below both the excellent 
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precision one can obtain in the primary cluster parameters (age, metallicity, distance, and 
reddening) and the range of related star cluster and stellar evolution problems that can be 
addressed. The goal of this paper is to present the Bayesian technique and demonstrate its 
internal precision. In subsequent papers we will derive improved WD and MSTO ages for 
clusters, with the long-term goal of intercalibrating WD and MSTO ages up to the ages of 
the oldest globular clusters. 

2. Baseline Stellar Cluster Model 

We chose a single set of stellar evolution ingredients to build and test the Bayesian 
approach. We use this model set to test the sensitivity of the derived WD and MSTO ages 
to the cluster parameters of [Fe/H], Ay, distance, age, number of cluster stars, and assumed 
photometric error. 

For our baseline Stellar Cluster Model we chose a Miller & Scalo (1979) initial mass 
function (IMF), main sequence and giant branch stellar evolution time scales of Girardi et 
al. (2000), the initial (main sequence) - final (white dwarf) mass relation of Weidemann 
(2000), WD cooling time scales of Wood (1992), and WD atmosphere colors of Bergeron 
et al. (1995). Using these ingredients we simulate star cluster color magnitude diagrams 
(CMDs), and using Bayesian techniques discussed below, we invert cluster CMDs to recover 
the probability distribution of the cluster parameters. 

When simulating a cluster, each star is randomly drawn from the IMF and, based on a 
user-specified binary star fraction, randomly assigned to be a single star or a binary with a 
companion also randomly drawn from the IMF. Note that although an IMF is required to 
simulate a cluster, the imphed age from either the MSTO technique or the WD technique 

is insensitive to the IMF. The IMF serves only to increase or decrease the population of 
stars of interest, e.g., MSTO stars or WDs. If there are insufficient stars, particularly if 
the cluster is young, then the few cluster stars coupled with the IMF can create a statistical 
uncertainty to locating the MSTO or perhaps even finding WDs. Binaries of nearly any mass 
ratio have a similar effect. WDs in binaries are generally not recognized and MSTO stars 
in such systems are found brighter and generally redder than the MSTO, and therefore they 
do not help define the MSTO. For these reasons and for simphcity in this study, we set the 
binary fraction to 0%, substantially lower than the typical value for open clusters of > 30%. 
For simplicity, we use only H-atmosphere (DA) WDs in our present simulations. While He- 
atmosphere (DB) WDs make up < 10% of field star WD samples (7% in Klcinman et al. 
2004), to date no DBs have been found in open clusters (Kalirai et al. 2005). A limitation of 
our cluster simulations is that stars with masses < O.25M0 are not included, thus producing 
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an unrealistic lower limit to the main sequence. Since the focus of this study is on stars that 
can become WDs, this simplification is merely one of presentation. 

Other stellar evolution (e.g., Yi et al. 2001; Baraffe et al. 1998; Siess, Dufour, & Forestini 
2000) and WD coohng (e.g., Benvenuto & Althaus 1999; Hansen 1999) models could have 
been used, and will be added to our code later, but for the present purposes, the above- 
mentioned, often-used models adequately cover parameter space and allow us to build and 
test the Bayesian machinery. 

After producing simulated CMDs we incorporate realistic photometric errors assuming 
reasonable cluster parameters, e.g., (m — My) = 12.5 and Ay = to 1, and assuming 
observations are obtained with the HST or similar imaging instrument able to observe to 
V = 27 with S/N = 15^. We use a conservative upper limit to the photometric precision of 
S/N=200, though we do not incorporate systematic calibration errors. Our Stellar Cluster 
Model limits are currently set by the Girardi et al. (2000) and Wood (1992) models, and these 
limits are 100 Myr to 4.5 Gyr and Z = 0.0004 to 0.030 ([Fe/H] « -1.676 to +0.198). This is 
adequate parameter space for significant age and metallicity exploration and to demonstrate 
the technique, though we clearly need to push the technique to greater ages. 

Our cluster simulations do not include mass segregation nor other dynamical processes, 
potentially important in open clusters, especially for the lowest mass stars; these typically 
have little effect on the measured WD mass fraction (von Hippel 1998; see also Hurley & 
Shara 2003 who find that the WD luminosity function and mass function are insensitive to 
dynamical effects at 0.5 to 1 half-mass radii). Simulated clusters specifically tuned to match 
real clusters using our Stellar Cluster Model have been presented in paper 1 (figs. 4-10). 
Here we do not attempt to match actual clusters, i.e., we do not tune distance, reddening, 
metallicity, cluster richness, and age, but rather we explore hypothetical clusters that cover 
the parameter space available to us. The CMDs for two such clusters are presented in Figures 
1 and 2 for log(age) = 9.0 and 9.5, respectively. The masses of a few WDs from across the 
cooling sequence are indicated in the second panels of both figures. Between 1 Gyr (Fig. 1) 
and 3.2 Gyr (Fig. 2) the WD terminus has evolved from My ~ 13 to My 14.5, and the 
simulated photometric errors have increased for the faintest WDs. 

In the next section, we outline the Bayesian technique that we will use in forthcoming 
studies to invert actual CMDs. In verifying the technique, rather than apply our Bayesian 
code to actual clusters with necessarily unknown parameters, we instead apply our code to 
simulated clusters. The analyses of simulated data sets test the degree to which an entirely 



^From experience. S/N=15 is required to obtain good morphological rejection of background galaxies at 
HST resolution (von Hippel & Gilmore 2000). 
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consistent set of stellar models, along with realistic photometric errors, yield the original 
input parameters. Our Bayesian analyses thus test the internal precision of our technique and 
its sensitivity to photometric errors, given the many non-linear aspects of stellar evolution. 
Since all stellar evolution models are imperfect, this approach provides a measure of internal 
precision only, not external accuracy. Our goal here is to build a modeling procedure with 
internal uncertainty < 5% in age, allowing us, when we subsequently analyze real clusters 
with high-quality data, to test for systematic problems in stellar models and ages of not 
much more than 5%. 



3. Bayesian Technique 

The goal of our Bayesian technique is to use information from the data and from our 
prior knowledge to obtain posterior distributions on the parameters of our model. Our prior 
knowledge is encoded in prior distributions on the model parameters. The model parameters 
include cluster parameters such as age and metallicity and an initial mass for each cluster 
star. These parameters are the inputs to our Stellar Cluster Model, which we use to derive 
predicted photometric magnitudes. The likelihood function then compares the predicted 
magnitudes with the observed (or simulated) data. 

Bayes Theorem relates the posterior distribution to the prior distribution and the like- 
hhood function, li M — (Mi, M2, . . . , Mjv) is a vector of initial masses of all stars in the 
cluster and = (T, [Fe/H], Ay, {m — My)) is a vector of cluster parameters, then we can 
treat our Stellar Cluster Model as a function G{M, 0) that maps every reasonable choice of 
{M, 0) to a resultant set of photometric magnitudes. To obtain the likelihood, we assume 
that the errors in our measurements arc independently distributed and Gaussian with known 
variance. Suppose there are stars in the cluster and we have observed them through n 
different filters. Then the observed data form an n X N matrix X with typical element Xij 
representing the magnitude in the ith filter of the jth star. By assumption, each magnitude 
is normally distributed: 

Xij N{fXij,a^j), (1) 

where fiij and afj are the mean and variance of the modeled photometry through filter i of 
star j. The means and variances also form nx N matrices, which we call /u and S. The full 
likelihood is then 

N / n 

p(xiM,s)=n n 

The variances S come from our knowledge of the precision of our observations. The means 
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fj, are the predicted photometric magnitudes that we obtain from the Stellar Cluster Model: 

l^ = G{M,@). (3) 

Thus, the hkelihood can be expressed in terms of the variables of our problem and the 
underlying Stellar Cluster Model: 

p(X|/x,S)=p(X|G,M,0,S). (4) 

Computationally, (2) is the most useful form of the likelihood because changing the under- 
lying Stellar Cluster Model leaves (2) unchanged. A different Stellar Cluster Model is just a 
different function, say H, such that 

lj = HiM,@) (5) 

or even 

lx^H{M,&). (6) 

for a different set of cluster parameters 0'. 

In Bayesian analysis, all model parameters require prior distributions. Wc have tried to 
select priors that are consistent with astronomers' knowledge of likely values for the various 
parameters. To reflect the fact that low mass stars are much more numerous than high mass 
stars and to be consistent with our Stellar Cluster Model where we used the Miller & Scalo 
(1979) IMF, we set the prior distribution on the logarithm of a star's mass proportional to 
the Gaussian distribution: 

p(log(M))o.exp( -('°«'^^' + ^-°^'\ (7) 

where the constants are from the fit derived by Miller & Scalo and the IMF is bounded at 
0.15 Mq and 100 Mq. For metaUicity, absorption, and distance modulus we use Gaussian 
priors in the common logarithm versions of these quantities ([Fe/H], Ay, {m — My)). We 
assume we have reasonable knowledge of the values and uncertainties of these parameters 
for a given cluster. This knowledge should come from outside information, not from the 
color-magnitude data we intend to analyze. Our prior on T, the base-10 logarithm of the 
cluster's age, is uniform between T = 8.0 and T = 9.7, and zero elsewhere. This is a power 
law prior on the age with exponent -1, which adequately reflects the observation that younger 
clusters are more common than older clusters. Note that priors from reliable, previously- 
derived cluster parameters are not required for our Bayesian approach, though they may 
help. The point is that priors encode any previously determined parameters, where they 
are available. In some cases constraining priors (e.g., small (j([Fe/H]) may turn out to be 
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required for precise results, in other cases, such as the ones studied here, constraining priors 
are unnecessary for precise results. 

Given the prior distributions and the likelihood, we obtain the posterior distributions of 
the parameters from Bayes theorem, which states that the posterior density p{9\y) on model 
parameters 9 given data y is 

(8) 

where p{y\6) is the likelihood and p{d) is the prior density on the model parameter 6. The 
denominator, is obtained by integrating the numerator over all possible values of 9 so 
that 

In our problem, it is impossible to compute the integral J p{y\9)p{9)d9 analytically. In- 
stead, we use Markov chain Monte Carlo (MCMC) to approximate the posterior distribution 
(Casella & George 1992; Chib & Greenberg 1995). The MCMC algorithm allows us to gen- 
erate a sample from the posterior distribution. We construct a Markov chain such that once 
it has converged, results of each iteration of the algorithm are distributed approximately 
according to the posterior distribution, and wc regard the history of the chain as a random 
sample from the posterior. We can thus obtain quantities of interest, such as sample means, 
without having to analytically compute the normalized posterior distribution. 

Our analysis relies on the Metropohs-Hastings algorithm (Chib & Greenberg 1995), 
which proceeds as follows: Suppose the current state at iteration t is 9^ = 9. Propose to 
move to some new state 9*. This proposal is generated with density q{9*\9). Compute the 
Metropolis-Hastings factor 

. \ p{9^\y)q{9\9^) J 
a = mm ' , , ' ' , , 1 (10) 
[pi9\y)qi9*\9) 'J ^ ^ 

and set 9*'^^ — 9* with probability a. Otherwise, set = 9. Our sample is the parameter 
sequence {9", . . . , 9^) where N is the total number of iterations and n is the number 
of iterations before the chain converges. We discard the first n — 1 iterations, which are 
referred to as the burn- in. Note the advantage of this method: since J p{y\9)p{9)d9 — 
J p{y\9*)p{9*)d9*, 

'p{y\9*)p{9*)q{9\9*) 



a — mm 



L P{y\9)p{9)q{9^\9) ' 

We can compute everything in (11) without calculating any intractable integrals. 



(11) 



The efficiency of the Metropolis-Hastings algorithm depends heavily on the choice of 
proposal distribution q. A common choice is a symmetric distribution centered at the cur- 
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rent value. This is the "random walk" Metropolis-Hastings sampler. This method has the 
advantages of simplicity and ease of implementation. However, the sampler can be inef- 
ficient if the distribution's width is inappropriate — the sampler might propose excessively 
small steps and take too long to traverse parameter space, or it might propose unreasonably 
large jumps and frequently reject steps. Another option is to choose a proposal distribution 
that approximates the posterior distribution. This kind of sampler is known as an "inde- 
pendence" sampler since q{0*\9) = q{6*), so that each proposed value is independent of the 
current state. The closer the proposal distribution approximates the target distribution, the 
higher the acceptance rate and (generally speaking) the more efficient the sampler. 



3.1. MCMC Sampling 

One of the chief problems in designing the MCMC sampler was overcoming the strong 
correlations between many of the variables. For instance, for a given position on the CMD, 
an increase in the age of the cluster will require a decrease in the mass of a WD, and vice 
versa. Since each parameter is sampled on individually in sequence, without removing the 
correlations the sampler can only take small steps in age or mass; if too large a step is taken, 
the proposed star's photometry will be too far from the observed position and the step will 
be rejected. In a similar way, metallicity and distance modulus are correlated with each 
main sequence star's mass, with each other, and with reddening. While from a theoretical 
standpoint removing these correlations is not required to obtain valid results, the number 
of iterations needed to be certain the entire posterior distribution is well sampled would 
necessitate far more computation time than is practical. 

Fortunately, over the ranges that our MCMC typically samples, these correlations are 
all nearly linear. To remove the WD age-mass correlation, we introduce a new parameter, 
U, and a constant, /3, defined by: 

Mk^P{n-T) + Uk, (12) 

where M^, Uk, and are the mass, dccorrelated mass parameter, and logarithm of the 
cluster age at the kth iteration, respectively, and T is the mean log cluster age. Then, rather 
than samphng on mass, we sample on U for each star. The MCMC algorithm then computes 
the mass at each iteration from the above equation. Figures 3 and 4 present the log(age) 
sampling history before and after correlation is removed, within the same MCMC run. In 
Figure 3, age values spanning ~ 100 iterations are correlated, meaning little new is learned 
about the posterior distribution within that correlation length. In Figure 4, the log(age) 
history is well-sampled and each iteration usefully samples the posterior distribution. The 
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new parameter, is then decorrelated from distance modulus and metallicity in a similar 
manner. Finally, the distance modulus and metallicity are decorrelated from one another 
and then from reddening. 

In order to improve the efficiency of our MCMC algorithm, we still need to address 
several sampling issues. For some parameters, the correlations become non-linear, often at 
their extreme values. For other parameters, the correlations consists of two or more separate, 
nearly linear, pieces with different slopes. For the brightest (youngest) WDs, the correlations 
between mass and age can be incredibly tight, and further work needs to be done for these 
objects to more precisely trace these correlations. 

The burn-in period for our MCMC runs consisted of a brief (5000 samples) period to 
settle close to the correct values and adjust step sizes. This was followed by two periods of 
5000 samples each to calculate the correlation between mass and age for WDs, two more to 
calculate the correlations between modulus and mass for main sequence stars, and between 
modulus and absorption, and two more to calculate the metallicity-main sequence mass and 
metallicity-absorption correlations. Finally, there is another 5000 sample period to adjust 
step sizes again. The whole burn-in, except for the initial settling-in period, is then repeated 
to more precisely determine the correlation factors, for a total burn-in period of 70,000 
samples. 

4. Demonstration and Discussion 

For the tests presented here we placed priors on cluster distance moduli, metallicities, 
and absorption values. The priors were normal distributions centered on the simulated Stellar 
Cluster Model parameters for (m — My), [Fe/H], and Ay, with the further requirement that 
A\ > 0. The Ay — runs are limiting cases, and for these we assumed cr(^v) = 0, i.e., 
there was no sampling on Ay. For the Ay = 0.1 and 0.3 cases, we assumed o"(Av) = 
0.1 and for the Ay = 1 case we assumed (j{Ay) = 0.3. For the other cluster parameters we 
assumed cr([Fe/H]) = 0.3 dex and a{m — My) = 0.2. All of these priors represent conservative 
uncertainties for well-observed, low-reddening clusters. We also placed a prior on the mass 
distribution for any given star with a form, as discussed above, based on the low-mass IMF. 

We simulated B, V, and / photometry for clusters for the range of parameters log(age) 
= 8.3, 8.7, 9.0, 9.3, and 9.5; [Fe/H] = -1.0, -0.15, 0.0, and +0.15; N (number of cluster 
stars fainter than the MSTO, including WDs) = 50, 100, 200, and 400; (m - My) = 12.5; 
and Ay — 0, 0.1, 0.3, and 1. Since our goal is to test the age sensitivity of the WDs, we 
removed all MSTO, subgiant, and giant branch stars from each simulation, so that what 
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remains are WDs and essentially unevolved main sequence stars. 

While we astronomers are most comfortable studying star clusters in the CMD, our 
Bayesian technique does not use the CMD, with its correlated errors between the x-axis and 
y-axis, but rather it uses an n-dimensional space, where n is the number of filters available 
and the units are magnitudes. In the numerical experiments we present here, n—3, as we 
use B, V, and / photometry. The input to the MCMC routine sees the CMD of Fig. 1, 
for example, in a form more akin to Figure 5, although offset by the simulated distance 
modulus. For presentation purposes, we reduced three dimensions to two by plotting either 
the B 01 I absolute magnitudes on the horizontal axis. One disadvantage of this plot is the 
large dynamic range in both axes. Still, the main CMD features can be discerned, e.g., WDs 
are clearly visible in the faintest corner of the plot. Reddening vectors for Ay = 1.0 are also 
shown, as is the effect of increasing distance modulus by 1.0 mag. Both distance and the 
reddening vectors are nearly parallel to the main sequence, especially the BV main sequence. 
Decreasing metallicity from [Fe/H] = 0.0 to —0.1 moves the main sequences in almost the 
opposite directions as the reddening vector. While there are some features in this diagram, 
and while the various distance, reddening, and metallicity vectors are not absolutely parallel 
and therefore not entirely degenerate, this diagram suppresses subtleties that primarily affect 
stellar color. 

Although we can simulate clusters younger than log(age) = 8.3, the MCMC technique 
requires sampling an age range, and for younger clusters this would often hit our (current) 
lower age limit of log(age) = 8.0. For the Ay cases, three simulated clusters were run 
for each unique set of cluster parameters. Any two clusters with identical parameters will 
yield different CMDs as both the IMF and the simulated photometric error distribution are 
sampled anew. After creating a cluster, we pass it to the MCMC routine with estimates 
of the mass of each star and estimates of the cluster parameters as starting points. (Our 
experiments show that as long as the MCMC algorithm converges, the results do not depend 
on the starting points. Starting points within a factor of ~ 2 in age or metallicity, for instance, 
are adequate for convergence.) Because the MCMC sampling is still often correlated, we 
sample for 10® iterations, reading out every tenth value for each stellar mass and for the 
cluster parameters. Many of our MCMC runs have a correlation length of < 10, and this 
produces uncorrelated parameter values. For those cases that still remain correlated, and 
guided by the rule-of-thumb that one typically wants 10^ uncorrelated iterations in order to 
adequately sample the posterior distribution, we find 10^ iterations works well for most of 
our simulated clusters. 

Figure 6 presents a well-sampled, typical history plot of cluster age for the cluster of 
Fig. 1. Figures 7 and 8 present the companion history plots for cluster [Fe/H] and (m — My). 
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There is a small amount of sticking in the sampling of these two variables at the beginning 
of the sequence, just after burn-in, and again near iteration 1.47 xlO^, but otherwise these 
history plots are well sampled. Since there is no Ay sampling in the Ay = case of Fig. 1, 
Figure 9 presents the Ay history plot for a cluster with the same parameters, except input 
Ay = 1. Histograms of these four types of history plots (Figure 10) are the estimates of the 
posterior probability distributions. In Figure 10 we present also the Ay = 0.3 and 1 cases. 
Figure 10 shows that the posterior distributions of log(age) and the other cluster parameters 
are close to normal, and furthermore that changing the absorption causes no strong bias 
in the results (more on absorption and bias below). From these posterior probabilities we 
can calculate statistics of interest, e.g., mean, median, a, percentiles, etc., and these are 
presented, below. Figure 11 presents the posterior probabilities of mass for four stars from 
the cluster simulation of Fig. 1. The first panel shows the mass posterior for a high mass 
WD, the second for a lower mass WD, the third for a main sequence star not far below 
the turn-off, and the final panel presents a low mass main sequence star. In all cases the 
mass distribution is nearly centered on the input mass-the mass value before its photometry 
was subject to random error-except in the case of the lowest mass star, where it differs by 
only 0.002 Mq. For the main sequence stars, the mass distribution is particularly narrow, 
showing that within the assumption of a specific model, precise photometry yields precise 
masses. The WD mass distributions are slightly broader because wc plot the zero-age main 
sequence (ZAMS) masses for these stars and a wider range of initial main sequence masses is 
converted into a narrow range of WD masses via the initial-final mass relation (Weidemann 
2000). 

In Figure 12, we check the differences between the mean log(age) values of the distri- 
butions and the input log(age) values compared to the standard deviations of the posterior 
log(age) distributions (cr). We are essentially asking what the deviation of each result is in 
units of its standard deviation. For the Ay = case, the distribution of errors is very similar 
to the overplotted normal distribution. Formally, the error distribution is also close to normal 
with average = 0.037, median = 0.055, standard deviation = 0.985, and skew = —0.191. This 
comparison is a sanity check on the self-consistency of our implementation of the Bayesian 
technique and whether the standard deviation statistic adequately captures the shapes of 
the posterior distributions. For the Ay — 0.1, 0.3, and 1 cases, our MCMC approach tends 
to be biased high in age by 0.56 a, 0.26 a, and 0.12 a, and for the higher absorption cases to 
have a pronounced non-normal distribution. These distributions are virtually identical if we 
plot median values instead of means of the posterior distributions. While we are still trying 
to understand some of the subtleties of the higher Ay cases, these offsets, which corresponds 
to 2.6%, 1.2%, and 0.6% systematic errors in age, are small enough that we set aside their 
resolution for now. Figure 12 demonstrates, particularly for low absorption values, that the 
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standard deviations in the posterior distributions are accurate assessments of the uncertain- 
ties due to photometric errors or any effects due to low number statistics, such as having 
very few WDs in young or sparse clusters. 

In Figure 13, we present the standard deviation log(age) uncertainty for each model 
with one or more WD versus the input log(age) of the cluster. The standard deviations are 
always small, typically < 0.04 dcx, corresponding to relative errors typically < 10%, for all 
ages tested. These errors do not depend significantly on the cluster age. In fact, the apparent 
slight dependence on age seen in Figure 13 is a combination of two other effects: there are 
fewer WDs in the youngest clusters and the coolest WDs in the oldest clusters are fainter, and 
therefore have higher photometric errors. Figure 14, which plots the same standard deviation 
uncertainties, now versus log(iVwD), shows the most important factor in this technique-the 
number of WDs per cluster. Although every WD contains age information, the quahty of 
that information is not the same for all WDs. Photometry for the coolest WDs in any cluster 
provides the most information (see below), yet photometric precision drops with decreasing 
WD luminosity. Therefore the age precision does not improve as the square root of the 
number of WDs, but somewhat more slowly. Nonetheless, even with 10 WDs the statistical 
(internal) error is almost always < 0.02 dex, or < 5%. Even with three WDs, 10% precision 
is usually achieved. 

Since the precision is so high via this technique, it is worth taking a small detour into 
the details of the age information locked up in each WD. Figure 15 presents the relationship 
between possible masses and possible ages for six of the nine WDs of Fig. 1 (log(age) = 
9.0). Two WDs are not presented because they are so close in mass to other WDs that 
they crowd the figure without presenting any new information. The faintest WD of Fig. 
1, has been dropped also, since its ZAMS mass (7.6 Mq) is beyond the 7.0 Mq limit of 
the Weidemann (2000) initial-final mass relation and its WD mass (1.131 Mq) is beyond 
the Wood (1992) 1.0 Mq cooling WD model hmit. Although this extrapolated WD has 
properties that are internally consistent enough for the MCMC runs, they are not pedagog- 
ically helpful in understanding the precision in the WD technique. Figure 15 shows that 
the age-mass relationship for the hotter, brighter WDs is highly correlated, which is the 
cause of the numerical difficulties with correlated sampling mentioned above. Is there useful 
age information in these hot, rapidly cooling WDs; information that is unexploited in the 
traditional approach, which uses the coolest WDs to derive a cluster's age? 

Figure 16 presents the age sensitivity of this technique^ by presenting the allowed age- 



^Note that the slight changes of slope in the mass-age relationship for the lowest mass WD is a numerical 
artifact cause by our use of linear interpolation among the Girardi et al. (2000) models. Besides being a 
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mass relationship for each of the WDs of Fig. 15, if each were the only WD in the otherwise 
identical cluster. Running the same cluster with only one WD yields the age constraints 
from an individual WD, while still relying on the cluster main sequence to constrain the 
combination of metallicity, distance, and reddening. Now it is clear that the higher mass 
WDs provide the tightest age constraints, eliminating significant age range allowed by the 
lowest mass, hottest WD, for example. Yet, in this case even the second hottest WD contain 
significant age information. Our technique can be pushed to the point where ages can be 
derived for clusters without observing the coolest WDs, and a companion paper will explore 
the sensitivity of that approach (Jeffery et al. 2006). The higher mass WDs have a much 
fiatter slope in the mass-age diagram since large changes in ZAMS mass do not appreciably 
change the contribution of precursor time scales, nor do they evolve at a much different rate 
as WDs, at least not in this age range. 

Figure 16 still begs the question: Why is there any age sensitivity when there is only 
one cluster WD? The short answer to this is that the WD region of the CMD is not highly 
degenerate. Though it may be possible to make a highly degenerate CMD with a combination 
of few cluster stars, high and uncertain reddening, uncertain metallicity, etc., generally this 
is not the case. There are also other constraints on the WD properties. A WD cannot have a 
mass higher than the upper limit for creating WDs (8 Mq in all our simulations, most likely 
7-9 Mq) and a WD cannot have a mass so low that stars with lower initial masses are still 
present on the main sequence. Changes in mass move a WD in the CMD along essentially 
the same vector as changes in age for hot WDs, where precursor ages are important, and 
in a perpendicular direction to age for cooler WDs, where precursor ages are unimportant. 
Figure 17 attempts to make this clear by plotting a small portion of the CMD of Fig. 1 
around the simulated WDs. The simulated WDs are presented as error bars, whereas the 
input values, before photometric scatter was added, are presented as filled circles. Here all 
cluster WDs, except the highest mass (7.6 Mq) simulated WD, arc plotted. The small '+' 
symbols connected by lines show the effect of holding mass constant while changing log(age) 
by ± 0.01 dex, or in the case of the two highest mass WDs, by changing log(age) by ± 
0.02 dex. Open squares show the effect of keeping log(age) = 9.0 while adjusting the ZAMS 
masses by ± 2%, or for the two highest mass WDs, by ± 5%. WD isochrones for log(age) 
= 8.9, 8.95, 9.0, 9.05, and 9.1 are overplotted. Ultimately, the Bayesian technique is so 
sensitive because minor changes in WD mass or cluster age of just a few percent move the 
expected location of any WD significantly in the CMD. Also, some types of photometric 



small effect and outside of the actual fit presented in Fig. 15, choppyness due to linear interpolation serves to 
occasionally slightly decrease the precision of our technique. Higher-order interpolations have not yet been 
necessary and would nearly double the run-time of our MCMC code. 
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error, e.g., the 2 a B — V color error of one of the WDs in the middle of the cooling 

sequence, cannot be matched by any realistic adjustment of cluster or stellar parameters, 
and thus this photometric error does not drive the fit for this WD. In Figures 15 and 16, for 
example, this WD is the object plotted third from the right. Its ZAMS mass (3.285 Mq) and 
age (log(age)=9.0) sit right in the center of the sampled values, so this color error had no 
meaningful effect. Errors in color could cause one to mistake a WD for a field star, however. 
The solution to this problem in real clusters with possibly contaminating field stars is better 
photometry or classification-level spectroscopy to confirm the WD. 

Besides deriving precise ages, the Bayesian technique also can derive precise values for 
the cluster parameters of metallicity, (m — My), and Ay. In all these cases the standard 
deviations in the posterior distributions are small, typically < 0.03 dex, < 0.02 mag, and < 
0.01 mag, respectively. All of these posterior uncertainties are an order of magnitude smaller 
than width of the prior distributions (0.3 dex, 0.2 mag, and 0.1-0.3 mag, respectively), 
demonstrating that high quality priors in these parameters are not generally needed, at 
least for low-to-moderate, single-valued absorption, and that the results arc insensitive to 
the exact assumed starting values for these parameters. The cluster photometry contains a 
wealth of information, and the Bayesian technique, along with an assumed model set, brings 
this out to high precision. For clusters with ten WDs, the age precision is typically better 
than 5%, easily meeting our needs for a precise statistical tool. 



5. Bayes Meets Star Clusters: Other Uses 

In our work to date, we have focused on cluster ages, particularly via the WD technique. 
Age via the MSTO technique will be next. Our MCMC code also derives values for the other 
major cluster parameters: metallicity, distance modulus, and line-of-sight absorption, along 
with the individual stellar property of mass. For our own purposes, we intend to upgrade 
our MCMC code to include binaries drawn from reahstic mass ratio distributions, field stars, 
DB (He atmosphere) WDs, and a wider range of standard stellar and WD evolution models, 
including ages up to globular cluster values. 

Here are a few further example uses for our Bayesian technique: 

1. In our effort to improve WD and MSTO ages we will systematically study main 
sequence and WD model parameters that affect ages, such as core convection prescriptions, 
non-standard elemental abundances, and diffusion in main sequence stars, as well as surface 
convection prescriptions and C/0 phase separation in cool WDs. 

2. We intend to study the sensitivity in the implied underlying parameters of simulated 
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and actual clusters to the initial-final mass relation as well as the upper mass limit for 
creating WDs. 

3. Our code derives mass posterior distributions for every object in the cluster. These 
mass estimates would be a good starting point for IMF studies, particularly since one can 
adjust the priors on mass to reflect different assumed IMFs and one can incorporate as input 
any stellar evolution model. By adjusting the prior on the IMF, one could see how many 
cluster stars are required before the resulting IMF is no longer sensitive to the prior. 

4. Once we add binaries to the MCMC code, we intend to study cluster binaries, 
their masses, and mass ratios. This would be a step forward from the typical approach of 
estimating cluster binary contributions by visually studying the distribution of stars above 
the single-star main sequence. Additional information, such as the probability of cluster 
membership from proper motion or radial velocity studies, could also be incorporated in the 
binary studies. 

We expect to make our code publicly available within a year, after it passes out of its 
development stage. 

6. Conclusions 

We have demonstrated a new Bayesian technique to invert color-magnitude diagrams 
to reveal the underlying cluster properties of age, distance, metallicity, and line-of-sight 
absorption, as well as individual stellar masses. Wc do not fit cluster fiducial sequences nor 
do wc create plots with many combinations of cluster parameters and then try to derive the 
best parameters via chi-by-cye. The Bayesian technique delivers not just parameters and 
error estimates, but entire posterior distributions. Posterior distributions for the parameters 
of interest are particularly valuable when they may be non-normal, as may occur with all 
the coupled, non-linear aspects of stellar evolution. Despite the potential for complex error 
distributions, we find posterior age distributions that are close to normal in log (age). Some 
other distributions, e.g., some mass distributions, are clearly non-normal. 

Within the confines of a given set of often-used models of stellar evolution, the initial- 
final mass relation, and WD cooling, and assuming photometric errors that one could rea- 
sonably achieve with HST, we find that our technique yields exceptional precision for even 
modest numbers of cluster stars. For clusters with 50 to 400 members and one to a few 
dozen WDs, we find typical internal errors of a"([Fe/H]) < 0.03 dex, a{'m — My) < 0.02 
mag, and a{Ay) < 0.01 mag. The parameter we are most concerned with, cluster WD age, 
has an internal error of typically only 0.04 dex (10%) for clusters with only three WDs and 
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almost always < 0.02 dex (< 5%) with ten WDs. All of these results have posterior distri- 
butions an order of magnitude narrower than the priors we applied, and therefore represent 
the actual information in cluster CMDs. Cluster photometry clearly contains a wealth of 
information, much of it coupled in a non-linear fashion, and the Bayesian technique, along 
with an assumed model set, brings this out to high precision. 

We thank the referee, Gordon Drukier, for helpful suggestions that improved our manuscript. 
This material is based upon work supported by the National Aeronautics and Space Admin- 
istration under Grant No. NAGS- 13070 issued through the Office of Space Science. We also 
gratefully acknowledge a grant to JS from the NSF-funded Vertical InteGration of Research 
and Education (VIGRE) program awarded by UT's Department of Mathematics. 
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Fig. 1. — (a) BV and (b) VI CMDs in the dereddened, absolute magnitude plane for a 
representative, log(age) = 9.0 cluster with [Fe/H] = 0.0, N = 100 main sequence plus WD 
stars, and photometric errors appropriate for (m — My) — 12.5, A\ — 0. Photometric errors 
in the WD region are similar in the BV and VI CMDs for these simulations, though there 
is a X-axis scale change. Representative ZAMS masses for the WDs are given in solar units. 

Fig. 2. — Same as Fig. 1, except for log(age) = 9.5. The oldest WDs are now fainter and 
have larger simulated photometric errors. 

Fig. 3. — Correlated age vs. ZAMS mass sampling for a WD during burn-in. 
Fig. 4. — The WD in Fig. 3, after age-mass decorrelation. 

Fig. 5. — The Bayesian code's version of the Fig. 1 CMDs. Filled squares represent the BV 
plane and open triangles represent the VI plane for solar metallicity stars. Cluster WDs 
are in the upper right, with the VI sequence shghtly above the BV sequence. The effects 
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of metallicity are shown by the open circles, which are the BV and VI main sequences, 
for [Fe/H] = —1. The offsets of distance and reddening are removed from this plot for 
presentation purposes. Reddening vectors for Ay = 1 for both planes are in the upper left, 
clS IS db vector showing the effect of increasing the cluster's distance modulus by 1 mag. 
The reddening vectors and distance offset vector are replotted near the main sequences to 
facilitate comparison. 

Fig. 6. — Typical history plot of cluster log(age), here for the 1 Gyr cluster of Fig. 1. In 
this case, the sampling was excellent. For clarity, only every 100th point is plotted, and only 
after the initial, 70,000 iteration burn-in period. 

Fig. 7. — Similar to Fig. 6, but for cluster [Fe/H]. This particular MCMC run had very mild 
correlation in [Fe/H]. 

Fig. 8. — Similar to Fig. 6, but for (m — My). 

Fig. 9. — Similar to Fig. 6, but for Ay in a cluster with Ay — 1. 

Fig. 10. — Histograms of the MCMC history plots, for (a) log(age) of Fig. 6, (b) [Fe/H] 
of Fig. 7, (c) (m — My) of Fig. 8, and (d), Ay of Fig. 9. For the first three panels the 
posterior distributions for Ay — 0, 0.3, and 1 are plotted and for the final panel the latter 
two distributions are plotted, since there was no sampling on Ay for the Ay — case. 
Every tenth iteration, which is the frequency with which these MCMC runs were read out, 
is incorporated into the histograms. 

Fig. 11. — Posterior probabilities of ZAMS masses for four stars from the cluster presented 
in Fig. 1. Panel (a) shows the mass posterior for a high mass WD, panel (6) for a lower mass 
WD, panel (c) for a main sequence star not far below the MSTO, and panel (d) for a low 
mass main sequence star. In all cases the means of the mass distributions are similar to the 
input masses, labeled with small vertical marks at the bottom of each panel. 

Fig. 12. — The differences between the mean ages of the posterior distributions and the input 
ages normalized by the standard deviations (a) of each posterior age distribution for the (a) 
Ay — 0, (6) Ay — 0.1, (c) Ay — 0.3, and (d) Ay — 1 cases. The distribution of errors is 
very similar to the overplotted normal distribution for the Ay — case and becomes subtly 
biased to +0.12 to +0.56 a for higher values of Ay. 

Fig. 13. — Derived standard deviations for each model with one or more WD versus the 
actual (input) ages of the clusters. The standard deviations are always small, typically 
< 0.04, corresponding to relative errors typically < 10%, for all ages tested. One result 
(among 319 runs) with (j(log age) = 0.105 at log(age) = 8.7 and A^wd = 1, is not plotted 
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for presentation purposes. 

Fig. 14. — Standard deviation uncertainties versus log(A'wD) for the same data as presented 
in Fig. 13. The precision of the age fit improves approximately as the fog of the number of 
WDs, which is an important factor under the observer's control. 

Fig. 15. — WD ZAMS masses versus cluster age for six of the nine WDs for the cluster in 
Fig. 1. The other three WDs are not plotted for clarity. 

Fig. 16. — WD ZAMS masses versus cluster age for six modified versions of the cluster 
presented in Fig. 1. In order to explore the mass-age correlations and see which WDs 
provide the greatest age constraints, nine clusters, each with only one WD from the original 
nine cluster WDs, were created. Again, only six of these mass-age relations are plotted for 
clarity. The lowest mass WDs have the tightest mass-age correlations, which creates greater 
MCMC sampling challenges. The higher mass WDs provide tighter age constraints. The 
kinks in the lowest mass relationship occur at boundaries of the main sequence (Girardi et 
al. 2000) tracks, and are numerical artifacts. 

Fig. 17. — WD regions of the Fig. 1 CMDs. The input WDs are plotted as filled circles and 
the scattered photometry are plotted as 1 a error bars. The highest mass WD is not plotted 
(see text). The symbols connected by lines show the effect of changing log(age) by ± 
0.01 dex, or in the case of the two highest mass WDs, by ± 0.02 dex. Open squares show 
the effect changing ZAMS masses by ± 2%, or for the two highest mass WDs, by ± 5%. 
WD isochrones for log(age) = 8.9, 8.95, 9.0, 9.05, and 9.1 are overplotted. The reddening 
vectors for Ay — 0.1 are also shown. 
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